# fit slugs models from lecture 18
slugs <- read.table( 'http://www.bio.ic.ac.uk/research/mjcraw/statcomp/data/slugsurvey.txt', header=TRUE)
 
glm.pois1 <- glm(slugs~1, data=slugs, family=poisson)
glm.pois2 <- glm(slugs~field, data=slugs, family=poisson)
 
library(MASS)
glm.NB1 <- glm.nb(slugs~1, data=slugs)
glm.NB2 <- glm.nb(slugs~field, data=slugs)
 
library(gamlss)
glm.NB3 <- gamlss(slugs ~ 1, sigma.formula=~field, data=slugs, family=NBI)
glm.NB4 <- gamlss(slugs ~ field, sigma.formula=~field, data=slugs, family=NBI)
 
LL.glm <- sapply(list(glm.pois1, glm.pois2, glm.NB1, glm.NB2, glm.NB3, glm.NB4), logLik)
k.glm <- sapply(list(coef(glm.pois1), coef(glm.pois2), c(coef(glm.NB1), glm.NB1$theta), c(coef(glm.NB2), glm.NB2$theta), c(coef(glm.NB3), coef(glm.NB3, what='sigma')), c(coef(glm.NB4), coef(glm.NB4, what='sigma'))), length) 
 
labels.glm <- c('glm.pois1', 'glm.pois2', 'glm.NB1', 'glm.NB2', 'glm.NB3', 'glm.NB4')
output.glm <- data.frame(model = labels.glm, LL = LL.glm, k = k.glm)
output.glm
 
# calculate AIC
AIC(glm.pois1, glm.pois2, glm.NB1, glm.NB2, glm.NB3, glm.NB4)
# add AIC to data frames
cbind(output.glm, AIC(glm.pois1, glm.pois2, glm.NB1, glm.NB2, glm.NB3, glm.NB4))
 
 
### Galapagos data 
gala<-read.table('ecol 563/galapagos.txt', header=T)
 
# normal model: richness as a function of Area
out.norm1 <- lm(Species~Area, data=gala)
# normal model: richness as a function of log Area
out.norm2 <- lm(Species~log(Area), data=gala)
# compare models
AIC(out.norm1, out.norm2)
 
# plot data
plot(Species~Area, data=gala)
abline(out.norm2, lty=2)
 
plot(Species~log(Area), data=gala)
abline(out.norm2, lty=2)
 
#Poisson model
out.pois1 <- glm(Species~log(Area),data=gala,family=poisson)
AIC(out.norm1, out.norm2, out.pois1)
 
# add Poisson model to graph
coef(out.pois1)
curve(exp(coef(out.pois1)[1]+coef(out.pois1)[2]*x), add=T, col=2, lty=2)
 
# identity link cannot work with log(Area) because
# some values are negative
out.pois2 <- glm(Species~log(Area), data=gala, family=poisson(link=identity))
# it does work with Area
out.pois2 <- glm(Species~Area, data=gala, family=poisson(link=identity))
AIC(out.norm1, out.norm2, out.pois1, out.pois2)
 
# negative binomial
library(MASS)
out.NB <- glm.nb(Species~log(Area), data=gala)
AIC(out.norm1, out.norm2, out.pois1, out.NB)
 
# add negative binomial to graph with Poisson
coef(out.NB)
curve(exp(coef(out.NB)[1]+coef(out.NB)[2]*x), add=T, col=4, lty=2)
 
# fit NB models with gamlss
library(gamlss)
 
# quadratic mean variance relationship--same as glm.nb
out.gamliss1 <- gamlss(Species~log(Area), data=gala, family=NBI)
AIC(out.NB, out.gamliss1)
# linear mean variance relationship
out.gamliss2 <- gamlss(Species~log(Area), data=gala, family=NBII)
AIC(out.NB, out.gamliss1, out.gamliss2)
 
# check mean variance relationship by grouping the response using values of Area
# quintiles of log(Area)
quantile(log(gala$Area), seq(0,1,.2))
 
# group observations into quantiles of Area
out.cut <- cut(log(gala$Area), quantile(log(gala$Area), seq(0,1,.2)), include.lowest=T)
 
# use area groupings to calculate mean and variance of richness
out.m <- tapply(gala$Species, out.cut, mean)
out.v <- tapply(gala$Species, out.cut, var)
 
# plot variance versus mean and superimpose quadratic fit
plot(out.v~out.m)
out.r <- lm(out.v~out.m+I(out.m^2))
coef(out.r)
curve(coef(out.r)[1]+coef(out.r)[2]*x+coef(out.r)[3]*x^2, add=T, col=2)